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Abstract 

We have carried out extensive isokinetic ah initio molecular dynamic simulations to investigate 
the finite temperature properties of the impurity doped cluster LigSn along with the host cluster 
Li7. The data obtained from about 20 temperatures and total simulation time of at least 3 ns is 
used to extract thermodynamical quantities like canonical specific heat. We observe a substantial 
charge transfer from all Li atoms to Sn which inturn weakens the Li-Li bonds in LigSn compared 
to the bonds in \a-j. This weakening of bonds changes the finite temperature behavior of LigSn 
significantly. Firstly, LigSn becomes liquid-like around 250 K, a much lower temperature than that 
of Liy (~ 425 K). Secondly, an additional quasirotational motion of lithium atoms appears at lower 
temperatures giving rise to a shoulder around 50 K in the specific heat curve of LigSn. The peak in 
the specific heat of Li7 is very broad and the specific heat does not show any premelting features. 
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I. INTRODUCTION 



Small clusters are known to behave differently than their bulk counterparts, mainly be- 
cause of their finite size. The finite size effects are reflected in the most of their properties, 
like equilibrium geometries, energetics, optical properties, ionization potential, polarizabili- 
ties, et c 112131415 . For example, equilibrium geometries of small as well as medium size clusters 
are considerably different than the bulk structures. Stability of a class of clusters depends 
on the filling up of geometric or electronic shell. The well known example of such a stability 
is that of alkali metal clusters displaying magic numbers. Further, the properties of homoge- 
neous clusters may be altered, some times quite significantly by adding an impurity. In the 
last few years some work has been reported on the impurity doped clusters^. These studies 
have brought out a number of interesting aspects like traping of an impurity, modifications 
in the equilibrium structures, changes in the bonding characteristics and enhancement or 
otherwise in the stability. A few studies have also been reported on the optical properties. 
Many of these properties have been shown to get influnced by the relative difference in the 
valence, ionic radii and electronegativity. 

A majority of the investigations on heterogeneous and homogeneous clusters pertain 
to the ground state properties. During the last few years a number of experimental^ 
and theoretica l 10 ! 11 ! 12 studies on finite temperature behavior of clusters have been reported. 
Finite temperature behavior of homogeneous clusters turns out to be very intriguing and 
has brought out many interesting features in contrast to the bulk behaviopi 3 -. For instance, 
clusters exhibit a broad peak in the specific heat which is sometimes accompanied by a 
premelting feature like a shoulder indicating possible isomarization. Recent experiment s£ 
and simulation a 11 ! 12 on tin clusters have shown that these clusters become liquidlike at 
higher temperatures than the bulk melting temperature (Tmbuik) contrary to the belief that 
solidlike to liquidlike transition in the finite size systems occurs at a temperature lower than 
Tmbuik because of the surface effects. A number of experiments on Na clusters by Haberland 
and co-workers^ have demonstrated that there is no systematic size dependence observed in 
the melting temperature of clusters between sizes 50 to 200. More interestingly, the peaks 
in the melting temperature as function of size do not occur at the geometric or electronic 
shell closing. 

In the present work, we investigate the finite temperature behavior of Li 6 Sn and Li 7 . 
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Since the impurity is known to change the geometry as well as bonding substantially in 
the host cluster it is also expected to change the finite temperature properties of the host 
cluster. The motivation for investigating these clusters comes from our previous study of Sn 
doped Li clusters^. This study has shown that, doping of lithium cluster with a tetravalent 
impurity such as Sn, induces a charge transfer from Li to Sn and as a result the nature of 
bonding changes from metallic like to dominantly ionic like. Preliminary investigation of 
the finite temperature behavior of LigSn has also indicated that this charge transfer induces 
dramatic effects at low temperatures. The interest in these systems also stems from the 
anomalous behavior seen in the bulk properties (e.g. electrical resistivity, density, stability 
function) of Li-Group IV alloys as a function of composition 15 . These alloys show peak in 
the resistivity around 20% Sn composition and this has been interpreted as atomic clustering 
around Sn at this composition whereas other alkali metal-Group IV alloys show peak in the 
resistivity around 50% composition of Group IV elements. It has also been observed that 
Li-Si alloys 16 show 15% volume contraction in the liquid state for Si concentration between 
15% to 50% whereas in Na-Sn alloys^ melting takes place via rotor phase in which alkali 
atoms are suppose to diffuse freely while the group IV polyanions perform only rotational 
or vibrational movement. 

As we shall demonstrate in this paper, the impurity induces a quasirotational motion of 
alkali atoms around Sn at low temperatures. In addition, in the liquid-like state the value 
of root mean square bond length fluctuations (5 rms ) of impurity doped cluster turns out to 
be considerably lower than that of Li7 indicating that volume expansion upon melting is 
much less in this cluster. Although it is difficult to establish a direct correlation between 
the properties of the bulk alloys and a small cluster like Li 6 Sn, we believe that the physics 
of impurity induced effects in both the systems essentially has the same origin arising out 
of the charge transfer which results into the strong ionic bond. 

It is now well established that for reliable ion dynamics one must reproduce the bonding 
properly. Description of ion dynamics critically depends on the correct microscopic descrip- 
tion of instantaneous state of electrons as system evolves in time. This is especially true 
when the nature of bonding may change during evolution of the trajectories at different 
temperatures. Therefore, for obtaining reliable ion dynamics we have chosen to use density 
functional theory(DFT) . Small size of the system has permitted us to perform ab initio cal- 
culations for substantially larger simulation time as compared to the previously published 
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ab initio result a 10 ^ 1 ! 12 . 

The paper is organized as follows. In Sec. |H] we describe briefly our computational and 
statistical approaches. Results for the equilibrium geometries, nature of bonding and finite 
temperature properties are given in Sec. 11111 and the conclusions are given in Sec. I1V1 

II. COMPUTATIONAL DETAILS 

All the calculations have been performed using the Kohn-Sham formulation of density- 
functional molecular dynamics (MD),— using the nonlocal norm- conserving pseudo- 
potentials of Bachelet et a/.— in Kleinman-Bylander form^fl, with the Ceperley-Alder— 
exchange-correlation functional. A cubic supercell of length 30 a.u. is used with energy 
cutoff of 15 Ry. for Li@Sn and 12 Ry. for LI7, which is found to provide sufficient con- 
vergence in total electronic energy(E tot ). To have a sufficiently precise evaluation of the 
Hellmann-Feynman forces on the ions, we ensure that the residual norm of each KS orbital, 
defined as I-HksV'i — £i.ipi\ 2 being the eigenvalue corresponding to eigenstate ipi of the 
KS Hamiltonian -£/ks)> was maintained at 10~ 12 a.u. After every ionic movement the total 
electronic energy is converged up to 10~ 5 a.u. for low temperatures (T< 200 K) and the 
convergence is 10~ 4 a.u. for higher temperatures. Such a high convergence is not warranted 
in the simulated annealing runs but it is very crucial for probing finite temperature be- 
havior of impurity doped clusters. The ground state and other equilibrium structures have 
been found by carrying out several steepest-descent runs starting from structures chosen 
periodically from a high-temperature simulation. 

The ionic phase space is sampled by isokinetic MD, 22 in which the kinetic energy is 
held constant using velocity scaling. For Li 7 we split the temperature range 100 K < T < 
1000 K into about 20 different temperatures, performing up to 90 ps of simulation for lower 
temperatures (i.e. T < 200) and 150-200 ps for higher temperatures. For LigSn we split 
the temperature range 10 K < T < 900 K into about 23 different temperatures, performing 
up to 200 ps of simulation for lower temperatures (i.e. T < 200) and 150 ps for higher 
temperatures. This results into total simulation time about 3 ns for LI7 and 5 ns for Li 6 Sn. 

The data obtained from these extensive simulations is used to extract different thermo- 
dynamic quantities such as entropy, specific heat, etc. The ionic specific heat is calculated 
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FIG. 1: Comparison of PED for 40 K and 175 K using different criteria for convergence in the 
total electronic energy (AE^Cf)- (a) denotes AEqqf = 0.00001 a.u. and (b) denotes Ai?scF = 
0.0001 a.u. 

from the fluctuations in the internal energy of the system and is given by, 

C(T) 



(1) 



where V is total electronic energy and (• • -) t denotes time averaged over the entire trajectory. 
A smooth specific heat can be obtained by using multiple histogram (MH) technique^. With 
multi histogram fit, which is basically a least square fit we extract density of states (Q(E)) 
or entropy {S(E) = ks In Q(E)). From entropy we can compute various thermodynamic 
quantities like canonical specific heat, C(T) = 6U(T)/dT, where U(T) = J E p(E,T) dE is 
the average total energy, and where the probability of observing an energy E at a temper- 
ature T is given by the Gibbs distribution p(E, T) = Q(E) exp(—E/k B T)/Z(T), with Z(T) 
the normalizing canonical partition function. 

Along with these thermodynamic quantities we also analyze the data using tradition 
parameters such as the £ rms , defined as 



2 ^ ((Rh)t-(Ru) 2 t ) l/2 

N(N - 1) f^ T 



(2) 



i>j (Ru)t 

where iV is the number of ions in the system, Ru is the distance between ions I and J, and 
(• ■ -)t denotes a time average over the entire trajectory. 
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We would like to make few pertinent comments on some technical issues which are specific 
to ab initio simulations. Since such DFT simulations are quite expensive, a judicious choice 
of simulation time step (At) , criteria for the convergence in the total electronic energy 
during Born-Oppenheimer MD (AE'scf) and total simulation time has to be made. Some 
of these issues may not be very crucial for geometry optimization, since there the details 
of ion dynamics are secondary. However for reliable extraction of thermodynamic data, the 
ion dynamics should be faithfully simulated using the accurate forces. Since the calculation 
of forces involves the total electronic energy (E tot ), the choice of AE SCF can turn out to be 
crucial. We demonstrate this point by showing potential energy distribution (PED) (shown 
in Fig. [TJ) calculated using two different AE SC f criteria, 0.0001 and 0.00001 a.u., for two 
temperatures 40 K and 175 K. Let us recall that PED(E,T) denotes accessible energy states 
between E and E + dE (where E is total electronic energy of the system). Fig. ^clearly 
shows that two different AE'scf leads to considerably different PEDs. The height of the 
peak in curve (a) {i.e. Ai^scF = 0.00001 a.u.) is about 25% higher than the curve (b) (i.e. 
AE'scf = 0.0001 a.u.) at 40 K. Further, it can be seen that the lower accuracy curve is seen 
to span a broader energy range and we have verified that this leads to the erroneous specific 
heat in this range of temperatures. However as the temperature rises lower criteria can be 
used without lose of accuracy. 

III. RESULTS AND DISCUSSIONS 

A. Finite temperature behavior of host cluster: Liy 

We begin our discussion by examining the equilibrium geometries of Li 7 (shown in Fig.|2J). 
The ground state structure is a pentagonal bipyramid with D 5 /j symmetry. The excited state 
is tricapped tetragon (C2 V symmetry) with 0.22 eV higher in energy compared to the ground 
state. These geometries are consistent with the earlier reported work^i As expected the 
charge distribution is completely delocalized indicating metallic like bonding in this cluster 
(not shown). 

Before discussing the thermodynamic quantities like specific heat it is instructive to ex- 
amine normalized potential energy distribution (PED) obtained via our molecular dynamic 
simulations, as a function of temperature. In Fig. 01 we show PED for different temperatures 
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AE = 0.0 eV 



E = 0.22 eV 



FIG. 2: The ground state and the first excited state geometry of Lij. AE represents difference in 
the total energies. 
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FIG. 3: Normalized potential energy distribution of Li7 for some relevant temperatures; viz 150 K, 
200 K, 240 K, 280 K, 320 K, 400 K and 500 K. It should be noted that PED for 400 K and 500 K 
are very fiat and almost overlapping 

ranging from 150 K to 500 K. It may be noted from Fig.|3]that the PED for low temperatures 
is very sharp and narrow indicating that the system spans a very restricted range of energies. 
Indeed the observed motion at these temperatures shows that atoms oscillate around their 
equilibrium positions. As temperature increases the pumped in energy is used to break some 
of the bonds and as a result cluster spans wider configuration space which reflects in the 
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FIG. 4: The canonical ionic specific heat of IA7. (a): Computed directly from fluctuations using 
Eq. ^with 75 ps for temperatures up to 200 K and 150 ps or more for higher temperatures, (b): 
The specific heat curve calculated using multihisto fit. 



broader PED. We note that around 400 K the PED becomes very flat which indicates that 
the cluster is spanning large number of energy states with equal probability. This rather 
broad PED around 400 K and above as compared to the PED for lower temperatures is 
suggestive of liquid-like nature of the cluster in this range of temperatures. 

Some insight can be gained by observing the detailed motion (the movie) of the cluster 
at few interesting temperature. We note that, around 360 K cluster starts revisiting the 
ground state via the first excited state and during this process, atoms from the pentagonal 
ring and atoms capping the ring get interchanged. This processes of interchanging of atoms, 
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FIG. 5: Root mean square bond length fluctuations of Li7 and LigSn. The slow monotonous rise 
in <5 rms for L17 below 320 K indicates solid-like behavior whereas linear nature above 450 K shows 
liquid-like behavior. For LigSn the Lindmann criteria gets satisfied at much lower temperature 
(around 175 K). We also note that for Li6Sn the value of <5 rms remains considerably constant for a 
wide range of temperature i.e. from 250 K to 700 K. 

becomes more frequent with increasing temperature. Indeed as we shall see this behavior 
leads to the peak in the specific heat around 425 K. In Fig. 01 we show ionic specific heat of 
Li7 calculated using Eq. directly (shown in Fig.HIa) as well as with the multiple histogram 
technique (shown in Fig.0jb). Evidently, both the methods yield similar specific heat except 
that specific heat by MH is much smoother. We note that the specific heat is very broad 
and does not change significantly after the peak around 425 K. Such a broad peak in the 
specific heat is a characteristic of small clusters in contrast to the sharp transition observed 
in the extended systems. 

These observations are also supported by the 5 rms which is shown in Fig. |S| As the name 
suggest <5 rms indicates fluctuations in the bond lengths averaged over all pairs and simulation 
time. According to the Lindmann criteria the clusters are considered to be in the liquid-like 
state when S Tms exceeds 0.1. It is interesting to note that <5 rms shows three distinct regions, 
a solid-like region below 320 K and a liquid-like region above 450 K. 320 K to 450 K is a 
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FIG. 6: The ground state and the first excited state geometry of LigSn. AE represents difference 
in the total energies. 

transition region and the value of <5 rms rises sharply in this region. If we chose to use the 
Lindmann criteria to estimate the "melting temperature" then this cluster "melts" between 
320 K to 360 K. 

To summarize, the specific heat curve indicates that above 425 K Li 7 is definitely liquid- 



like and the transition is characterized by rather broad peak in t 



re specific heat. These 



features are consistent with the specific heat of Li 8 reported in Ref. 
namely Li 7 and Lig do not show any premelting features. 



Both these clusters 



B. Finite temperature behavior of impurity doped cluster: LigSn 

In this section we present our results on the impurity doped cluster, namely LigSn and 
contrast the behavior with that of host cluster Li 7 . 

Before studying these finite temperature effects let us examine the equilibrium structures 
of Li 6 Sn which are shown in Fig. [HI Firstly, we note that both the geometries are considerably 
different than that of Li 7 . Secondly, in the ground state, Sn has six fold coordination whereas 
the first excited state which is 0.16 eV above the ground state, has five fold coordination 
with Li atoms. Thus the first excited state has one less Li-Sn bond. Further, these two 
structures can be obtained by capping the most stable I^Sn tetrahedron^ in two different 
ways. 

Now we turn to the bonding features of the system. Our earlier work 14 has shown that a 
tetravalent impurity like Sn in lithium clusters changes the delocalized metallic like bonding 
to dominantly ionic like. In Fig. [7] we show the difference charge density of Li 6 Sn, defined as 
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(a) (b) 




FIG. 7: The isosurfaces for difference charge-density of LigSn. (a): charge gained region, (b): 
charge depletion region. Sn is shown by the black sphere at the center whereas white spheres 
indicate Li atoms. 




abed 

FIG. 8: Quasirotational motion of Li atoms around Sn; (a): the ground state geometry (b and 
c): two Li-Li bonds are broken and the Li atoms get separated into two bunches which form a 
quasirotational motion around Sn. (d): The resultant geometry. We note that (a) and (d) denotes 
same geometries with different orientation. 

the difference between the selfconsistent charge density and the superimposed atomic charge 
density. Quite clearly, the charge gained region is mainly around Sn whereas the charge 
depletion region is near all Li atoms. This clearly shows that there is a charge transfer from 
all Li atoms to Sn, filling in the Sn centered p orbitals^S. Typically the depleted charge is 
of the order of 2/3 electron per Li atom. As a consequence of the significant charge transfer 
from Li to Sn, considerable depletion of charge in the region of Li-Li bond is observed which 
results into weakening of Li-Li bond. This weakening of Li-Li bonds show dramatic effects 
on the finite temperature behavior of the cluster. 

The consequences of weakening of Li-Li bonds are observed at very low temperatures. 
Around 30 K two Li-Li bonds break and Li atoms get separated into two bunches (see Fig. |H] 
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FIG. 9: Normalized Potential Energy Distribution at 150 K and 400 K of LigSn and Li7. 

b and c ). These bunches perform a quasirotational motion around Sn as shown in Fig. [HI We 
note from Fig. |H] that the initial (Fig. [H] a) and final (Fig. [H] d) geometries are same but have 
different orientations. Since the motion is observed at very low temperatures the potential 
barrier joining these two (identical) structures must be very low. This quasirotational motion 
seen at such a low temperature shows up in the PED of Li 6 Sn. In Fig.|§]we compare PED of 
LigSn with that of LI7. Evidently, at 150 K the PED of LigSn is significantly different than 
that of Liy. A rather narrow and sharp distribution for Li7 indicates that it is showing clear 
solid-like behavior exhibiting small oscillations at this temperature, whereas the significantly 
broader distribution of LigSn shows that other modes of excitations are also available to the 
system at this temperature. In fact the observed motion at this temperature shows breaking 
of most of the Li-Li bonds. On the other hand at 400 K the PED of both clusters have 
similar shapes. 

Turning to the ionic specific heat of LigSn (shown in Fig. 110)1 we note that the specific heat 
curve is significantly different than that of Li7. Unlike Li 7 , it shows a shoulder around 50 K— 
and the main peak is observed at substantially lower temperature (i.e. 250 K ) compared 
to that of Li7 (425 K). The shoulder in the specific heat curve of LigSn corresponds to the 
quasirotational motion described earlier. As temperature rises further this quasirotational 
motion becomes more rigorous and eventually around 200 K the bunches of Li atoms lose 
their identity and exhibit diffusive motion around Sn. This leads to the peak around 250 K 
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FIG. 10: The canonical ionic specific heat of LigSn. The shoulder around 50 K is because of the 
quasirotational motion seen at that temperature whereas the peak is around 250 K, a substantially 
low temperature compared to that of Liy. 

in the specific heat. 

These observations are also consistent with the <5 rms which is shown in Fig. |3] along with 
the 5 rma for Liz. For LigSn the value of <5 rms rises in two distinct steps. The first step is 
around 30 K. We recall that quasirotational motion starts around this temperature and 
another step is around 200 K where the bunches of Li atoms lose their identity and diffuse 
around Sn. 

It is also observed that the value of <5 rms remains nearly constant between 250 to 700 K. 
Although this effect is small it is noteworthy. In this temperature range the system is liquid- 
like but instead of random diffusive motion of all atoms in the cluster we observe that Li 
atoms perform constrained diffusive motion with Sn at the center. The averaged fluctuations 
in Li-Sn bond in this range of temperatures is found to be within 15%. Thse observations 
indicate that the main peak around 250 K is associated with the breaking of weak Li-Li 
bonds. We also note that in Li 6 Sn though <5 rms satisfies the Lindmann criteria at very low 
temperature compared to that of Li 7 , its saturated value is 0.17 which is considerably lower 
compared to the host Li 7 (5 rms = 0.25) cluster. This indicates that in the liquidstate the 
volume expansion of the impurity system is much less than in the case of host cluster. 
This is consistent with the existence of stronger Li-Sn bond in the system. At still higher 
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temperatures one of the ionic bonds breaks and we observe the first excited state (shown in 
Fig. El) in which Sn has one less coordination number than the ground state. Thus the route 
for evaporation proceeds via the first excited state. 

Finally we note that, Li 6 Sn is an over lithiated cluster, i.e. the number of Li atoms in the 
cluster are more than those required to complete the p shell of Sn atom. It will be interesting 
to investigate the finite temperature behavior of a cluster which satisfy the octate rule i.e. 
like Li4 X Sn x . In these clusters though we expect that because of complete charge transfer 
from all Li atoms to Sn, the bonds between Li atoms will be extremely weak, the ionic bond 
between Li-Sn will be more strong than those in the overlithiated cluster like LigSn. 

IV. CONCLUSIONS 

We have presented an ab initio isokinetic MD simulations on LigSn and Li7 to demonstrate 
the dramatic effects induced by the impurity in the host cluster. The key point which emerges 
out of these calculations is that the essential physics governing these effects comes from the 
charge transfer from Li to Sn because of the large difference in their electronegativities. In 
LigSn, we observe a quasirotational motion at very low temperature which gives rise to a 
shoulder in the specific heat curve around 50 K. The main peak is observed at much lower 
temperature, around 250 K, compared to the host cluster (« 425 K) whereas Li7 shows a 
broad peak in the specific heat with no premelting features. The motion in the liquid state 
of Li 6 Sn is essentially a constrained motion of Li atoms around Sn center. 
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